Introduction

 

Pufferfish, Takifugu rubripes is one of marine commercial fishes in Asia (Hamasaki et al. 2017; Kong et al. 2019). This species has high nutritional and economic value. In China, its production is approximately 70% of the total global production per year, and most are exported to Japan and Korea (Guo et al. 2017). However, various diseases caused by Edwardsiella tarda, Flexibacter maritimus, Myxococcus sp. and red sea bream iridovirus, have hindered its production and caused economic losses (Cui et al. 2014a; Ma et al. 2014). Therefore, it is urgent and significant to understand molecular immune-related mechanisms of T. rubripes.

RNA-Seq sequencing has developed rapidly over the past decade and has been widely used for transcriptome expression profiling (Zhou et al. 2016a, 2019), regulatory network construction (Zhao et al. 2016; Wang et al. 2018), noncoding RNA identification (Cui et al. 2017, 2020; Jiang et al. 2018), genetic diversity determination (Cui et al. 2014a; Zhou et al. 2016b) and other applications. RNA-Seq allows the identification of many immune-related genes and elucidation of mechanisms of immune evasion in fishes. For example, after infection with Aeromonas hydrophila, 2,900 differentially expressed genes (DEGs) are identified from common carp spleen, and these genes are mainly involved in immune response (Jiang et al. 2016). During channel catfish-Flavobacterium columnare interaction, the RNA-Seq results have shown that fish gill is involved in the immune response to the pathogen (Peatman et al. 2013; Sudhagar et al. 2018). Meanwhile, some defense-related genes and pathways have been identified form Lateolabrax japonicus and Cynoglossus semilaevis challenged with Vibrio anguillarum (Zhang et al. 2015; Zhao et al. 2016), Larimichthys crocea, Megalobrama amblycephala and Ctenopharyngodon idella challenged with A. hydrophila (Mu et al. 2010; Tran et al. 2015; Yang et al. 2016; Song et al. 2017), tilapia challenged with Streptococcus agalactiae and S. iniae (Zhu et al. 2017; Wang et al. 2016) and zebrafish challenged with Salmonella typhimurium, E. tarda and Mycobacterium marinum (Stockhammer et al. 2010; Ordas et al. 2011; Yang et al. 2012; Chen et al. 2017) by RNA-Seq.

In fishes, immune organs include liver, spleen, kidneys and gills. In our previous study, many immune-related genes and pathways have been found in the gills of T. rubripes (Cui et al. 2014b), and it is also found that these immune-related genes contain a number of single-nucleotide polymorphisms (SNPs) (Cui et al. 2014c). In addition, fish livers, spleens and kidneys, as important immune organs, play important roles in the response of fish against various pathogens. Chemokines are key immune regulators, acting as bridges between innate and adaptive immunity. Some catfish β and α chemokines are responsive to E. ictaluri infection in the lever (Fu et al. 2017a, b). Meanwhile, fish toll-like receptors, which are important immune-related proteins, can recognize pathogen-associated molecular patterns. During A. hydrophila challenge, the expression levels of genes encoded toll-like receptors are significantly changed in the kidneys and spleen of common carp, suggesting that the toll-like receptor genes affect in the host immune response to A. hydrophila infection (Gong et al. 2017).

Previous studies have shown that immune responses are different between males and females, which likely contributes to differences in the response of the two sexes to pathogens (Pennell et al. 2012; Klein et al. 2015). Females exhibit lower burdens of pathogen infection (Pennell et al. 2012). In previous studies, it has been found that a number of immune-related genes are differentially expressed between males and females. For example, the expression levels of some cytokines genes, including TNF-α, IL-1β, IL-6 and CXCL10, are more up-regulated in male mice than in females (Kahlke et al. 2000; Asai et al. 2001; Moxley et al. 2002; Marriott et al. 2006; Klein et al. 2015). However, the studies on the differences between the immune responses of male and female fishes are limited.

In this study, RNA-Seq is performed to investigate the differences in genome-wide gene expression in important organs, namely, livers, spleens and kidneys, between male and female T. rubripes. Sets of DEGs were identified. The important pathways and key immune-related genes are screened by Gene Ontology (GO), KEGG pathway enrichment and protein-protein interaction (PPI) network analysis. These results will not only provide insight into the molecular immune-related mechanisms of T. rubripes but also help future molecular breeding.

 

Materials and Methods

 

Ethics statement

 

All the fish experiments were conducted in accordance with the Guidelines for Experimental Animals established by the Ministry of Science and Technology. This study was approved by the Animal Care and Use Committee of the College of Fisheries and Life Science, Dalian Ocean University. All surgery was performed under sodium pentobarbital anesthesia, and all efforts were made to minimize suffering.

 

Sample collection

 

A total of 20 T. rubripes (10 males and 10 females) were sampled from Dalian Tianzheng Industrial Co., Ltd. (Dalian, China). Each individual was 15–18 cm long. T. rubripes were transferred to experimental tanks. These individuals were kept in sea water at 20 for 7 days to acclimate them to the experimental conditions before the study. The livers, spleens and kidneys of these fish were collected. The same tissue types from male or female fish were pooled. The tissues were placed in RNAlater (Ambion), stored at room temperature for 24 h, and then moved to -80 for storage until RNA isolation.

 

RNA extraction, library construction and sequencing

 

Total RNA of the livers, spleens and kidneys was extracted using RNAiso Plus (TaKaRa) according to the manufacturer’s instructions. Total RNA was sent for next generation-sequencing to Biomarker Biotechnology Corporation (Beijing, China). The library was sequenced on an Illumina HiSeq 2000 with 101-bp paired-end reads.

 

Read mapping

 

Based on the method of Kong et al. (2018), the low-quality reads were removed. The clean reads mapped to the T. rubripes genome (https://www.ncbi.nlm.nih.gov/genome/63?genome_assembly_id=22739) by TopHat.

 

Identification of new genes

 

To improve the annotation of the T. rubripes genome, we used the RNA-Seq data to identify new genes. First, Cufflinks was used to assemble all the mapped reads and map them to the genome. Then, sequences encoding short peptide chains (less than 50 amino acid residues) and those containing single exons were removed. Lastly, the new genes identified above were used as query sequences to search the Pfam, NR, eggnog, KOG, Swiss-Prot, GO, KEGG databases, and COGusing BLASTX.

 

Analysis of DEGs

 

The expression levels of the genes were measured as read counts normalized to the respective lengths of the genes using Cufflinks. DEGs of nine groups, namely, ♂liver vs. ♂spleen, ♂liver vs. ♂kidney, ♂spleen vs. ♂kidney, ♀liver vs. ♀spleen, ♀liver vs. ♀kidney, ♀spleen vs. ♀kidney, ♂liver vs. ♀liver, ♂spleen vs. ♀spleen and ♂kidney vs. ♀kidney, were identified using DEGseq with a false discovery rate (FDR) < 0.01 and fold change > 2 or < -2.

 

Gene ontology (GO) and KEGG enrichment analyses

 

We used the GOseq R packages to conduct GO enrichment analysis of the DEGs (Young et al. 2010), and KOBAS software was used to perform the KEGG pathway enrichment analysis (Mao et al. 2005; Kanehisa et al. 2008).

 

Protein-Protein interaction (PPI) analysis

 

The sequences of the DEGs were compared by BLASTX against the genome of a related species (information about the PPIs of which can be found in the STRING database: http://string-db.org/) to obtain the predicted PPIs associated with these DEGs. Then, the PPIs associated with these DEGs were visualized using Cytoscape (Shannon et al. 2003).

 

Results

 

RNA-Seq data from the livers, spleens and kidneys of T. rubripes

 

RNA-Seq was performed on six RNA samples (♂liver, ♂spleen, ♂kidney, ♀liver, ♀spleen and ♀kidney) from T. rubripes. A total of 274,074,082 clean reads were generated (Table 1). We obtained clean reads for approximately 41 Gb in total. The number of clean sequences from each library ranged from 41,722,196 to 52,936,660. The Q30 value of each library was more than 92%. The accession number of all clean data in the NCBI were SRR6925615, SRR6925526, SRR6925454, SRR6925406, SRR6925404 and SRR6925377.

 

 

Fig. 1: Unique reads mapped to various genomic regions

 

 

 

Fig. 2: Scaling normalization box plot of the gene mapping counts in the six samples. The similarity in transcript expression and variability of the six corresponding transcriptome libraries indicates that these libraries were comparable for the identification of DEGs at the transcriptome level. The y-axis represents the logarithm of the FPKM value. The x-axis represents the samples

 

Mapping to T. rubripes genome and identification of new genes

 

TopHat was used to align all the clean reads to T. rubripes genome (Trapnell et al. 2009). The mapping ratio ranged from 81 to 84% (Table 1). The number of unique mapped reads from ♂liver, ♂spleen, ♂kidney, ♀liver, ♀spleen and ♀kidney samples were 42,013,959, 34.293,686, 35,711,201, 33,536,438, 36,814,376 and 35,030,165, respectively. Sets of unique mapped reads were mapped to exons (88.09%/84.76%/86.43%/88.68%/85.39%/85.39%, ♂liver/♂spleen/♂kidney/♀liver/♀spleen/♀kidney), introns (6.17%/8.48%/7.03%/5.9%/7.97%/7.8%, ♂liver/♂spleen/♂kidney/♀liver/♀spleen/♀kidney) and intergenic regions (5.74%/6.76%/6.54%/5.42%/6.64%/6.81%, ♂liver/♂spleen/♂kidney/♀liver/♀spleen/♀kidney) (Fig. 1).

In addition, after analysis of the assembly and mapping data, we identified 1,443 new genes which were annotated by a BLASTX search against eight databases, namely, NR, Swiss-Prot, GO, COG, KOG, Pfam, eggNOG and KEGG. A total of 1,041 new genes had significant hits in at least one database.

 

DEGs analysis

 

The FPKM value of each gene was used to represent the normalized expression value. Scaling normalization of the count data indicated equal variation in the data for the transcriptomes of the six samples, as indicated in Fig. 2. A total of 3,673, 3,919, 2,912, 3,912, 3,840, 2,971, 369, 323 and 263 genes showed significantly differential expression in nine groups, namely, ♂liver vs. ♂spleen, ♂liver vs. ♂kidney, ♂spleen vs. ♂kidney, ♀liver vs. Table 1: Statistical data of RNA-Seq reads for six samples

 

Sample

Clean reads

Clean bases

Q30 (%)

Mapped Reads

Unique Mapped Reads

♂liver

52,936,660

7,920,730,888

92.55

44,708,721 (84.46%)

42,013,762 (79.37%)

♂spleen

43,129,860

6,439,115,608

92.58

35,581,354 (82.50%)

34,293,686 (79.51%)

♂kidney

45,024,504

6,723,621,318

92.53

36,914,930 (81.99%)

35,711,201 (79.32%)

♀liver

41,722,196

6,248,786,778

92.65

35,299,378 (84.61%)

33,536,438 (80.38%)

♀spleen

46,960,728

7,033,459,874

92.27

38,094,434 (81.12%)

36,814,376 (78.39%)

♀kidney

44,300,134

6,626,705,206

92.44

36,099,012 (81.49%)

35,030,165 (79.07%)

 

Table 2: Number of differentially expressed genes in immune organs

 

DEG set

DEG number

Up-regulated

Down-regulated

♂liver vs. ♂spleen

3,673

1,067

2,606

♂liver vs. ♂kidney

3,919

1,691

2,228

♂spleen vs. ♂kidney

2,912

1,855

1,057

♀liver vs. ♀spleen

3,912

1,417

2,495

♀liver vs. ♀kidney

3,840

1,830

2,010

♀spleen vs. ♀kidney

2,971

2,037

934

♂liver vs. ♀liver

369

172

197

♂spleen vs. ♀spleen

323

137

186

♂kidney vs. ♀kidney

263

175

88

 

 

 

Fig. 3:MA (log ratios–mean average) plots showing gene expression. The y-axis represents the log fold change, and the x-axis represents the log transcript counts. M is the log ratio of the two dyes used in the hybridization experiment, and A is the average of the log intensities. Red, green and black points represent upregulated, down-regulated and normally expressed genes, respectively

♀spleen, ♀liver vs. ♀kidney, ♀spleen vs. ♀kidney, ♂liver vs. ♀liver, ♂spleen vs. ♀spleen and ♂kidney vs. ♀kidney (Table 2), in which 3,555, 3,840, 2,844, 3,815, 3,754, 2,911, 352, 307 and 252 DEGs, respectively, were annotated by a BLASTX search against eight databases. The Fig. 3 showed that MA plots were drafted using “eps” format files. Compared to ♂spleen and ♂kidney, 1,067 and 1,691 DEGs, respectively, were upregulated and 2,606 and 2,228 DEGS, respectively, were downregulated in ♂liver. Between ♂spleen and ♂kidney, 1,855 DEGs were upregulated in ♂spleen and 1,057 DEGs were upregulated in ♂kidney. Meanwhile, compared with ♀spleen and ♀kidney, 1,417 and 1,830 DEGs were upregulated, respectively, and 2,495 and 2,010 DEGs were downregulated, respectively, in ♀liver. In ♀spleen vs. ♀kidney, 2,037 DEGs were upregulated in ♀spleen and 934 were upregulated in ♀kidney.

In addition, compared to female T. rubripes, 172, 137 and 175 DEGs were upregulated in the livers, spleens and kidneys of male T. rubripes. A Venn diagram of the DEGs showed that a majority of these genes were not shared by these nine groups. Only 553 DEGs were present in ♂liver vs. ♂spleen, ♂liver vs. ♂kidney and ♂spleen vs. ♂kidney, and 621 DEGs were present in ♀liver vs. ♀spleen, ♀liver vs. ♀kidney and ♀spleen vs. ♀kidney (Fig. 4). In addition, the DEGs were not shared between males and females; 323, 271 and 210 DEGs were present in only ♂liver vs. ♀liver, ♂spleen vs. ♀spleen and ♂kidney vs. ♀kidney, respectively (Fig. 4); 547 were found in only ♂liver vs. ♂spleen; and 786 DEGs were found in only ♀liver vs. ♀spleen (Fig. 4). Meanwhile, 808 and 710 DEGs were found in only ♂liver vs. ♂kidney and ♂spleen vs. ♂kidney, respectively, and 729 and 739 DEGs were found in only ♀liver vs. ♀kidney and ♀spleen vs. ♀kidney, respectively (Fig. 4).

 

GO enrichment analysis

 

GO enrichment analysis was performed, and the most significant enriched GO terms under Biological Process (BP) and Cellular Component (CC) were “G-protein coupled receptor signaling pathway (GO:0007186)” and “integral component of membrane (GO:0016021)” in ♂liver vs. ♂spleen, ♂liver vs. ♂kidney and ♂spleen vs. ♂kidney. On the other hand, the most significant enriched GO term under Molecular Function (MF) was “G-protein coupled receptor activity (GO:0004930)” in ♂liver vs. ♂spleen and ♂spleen vs. ♂kidney and “olfactory receptor activity (GO:0004984)” in ♂liver vs. ♂kidney. The enriched BP GO terms “one-carbon metabolic process (GO:0006730)” and “alpha-amino acid catabolic process (GO:1901606)” were present in only ♂liver vs. ♂spleen, and “regulation of cyclin-dependent protein serine/threonine kinase activity (GO:0000079)” and “killing of cells of other organism (GO:0031640)” were present in only ♂spleen vs. ♂kidney. The CC GO terms “postsynaptic membrane (GO:0045211)” and “connexin complex (GO:0005922)” were significantly enriched in ♂liver vs. ♂spleen and ♂spleen vs. ♂kidney, respectively. The enriched MF GO term “electron carrier activity (GO:0009055)” was present in only ♂liver vs. ♂spleen, “cation-transporting ATPase activity (GO:0019829)” in only ♂liver vs. ♂kidney, and “structural constituent of cytoskeleton (GO:0005200)” and “voltage-gated potassium channel activity GO:0005249” in only ♂spleen vs. ♂kidney.

Meanwhile, in female samples, the most significant enrichment GO term of BP and CC was also “G-protein coupled receptor signaling pathway (GO:0007186)” and “integral component of membrane (GO:0016021)”. While the most significant enrichment GO term of MF “G-protein coupled receptor activity (GO:0004930)” was in ♀liver vs. ♀kidney and ♀spleen vs. ♀kidney, and olfactory receptor activity (GO:0004984) was in ♀liver vs. ♀spleen. The enrichment GO terms of BP “alpha-amino acid catabolic process (GO:1901606)” and “transmembrane transport (GO:0055085)” were enriched in only ♀liver vs. ♀spleen and ♀spleen vs. ♀kidney, respectively; the GO terms of CC “proton-transporting V-type ATPase complex (GO:0033176)” and “intermediate filament (GO:0005882)” were only present in ♀liver vs. ♀kidney and ♀spleen vs. ♀kidney, respectively; The enrichment GO terms of MF “electron carrier activity (GO:0009055)” and “endopeptidase inhibitor activity (GO:0004866)” were only present in ♀liver vs. ♀spleen. “NADP binding (GO:0050661)” in only ♀liver vs. ♀kidney and “structural constituent of cytoskeleton (GO:0005200)” and “voltage-gated potassium channel activity GO:0005249” in only ♀spleen vs. ♀kidney.

The enriched GO terms were also compared between male and female T. rubripes. The enriched BP GO terms were not differentially expressed in ♂liver vs. ♀liver, ♂spleen vs. ♀spleen and ♂kidney vs. ♀kidney. The enriched CC GO term “myosin complex (GO:0016459)”and MF GO term “voltage-gated potassium channel activity (GO:0005249)” were present in only ♂spleen vs. ♀spleen.

 

KEGG pathway enrichment analysis

 

The top 20 pathways identified from each group by KEGG enrichment analysis are shown in Fig. 5. In ♂liver vs. ♂spleen, ♂liver vs. ♂kidney and ♂spleen vs. ♂kidney, the q-values of 20, 18 and 7 pathways, respectively, were lower than 0.05. The pathway categories “cytokine-cytokine receptor interaction” and “glycine, serine and threonine metabolism” were identified as being significantly enriched in ♂liver vs. ♂spleen, ♂liver vs. ♂kidney and ♂spleen vs. ♂kidney. “Fatty acid degradation”, “fatty acid metabolism”, “glycerolipid metabolism” and “tryptophan metabolism” were identified as significantly enriched categories in only ♂liver vs. ♂spleen; “one-carbon pool by folate” and “renin-angiotensin system” were significantly enriched in only ♂liver vs. ♂kidney; “carbon metabolism”, “ECM-receptor interaction”, “pentose phosphate pathway”, “neuroactive ligand-receptor interaction”, “pentose phosphate pathway” and “taurine and hypotaurine metabolism” were significantly enriched in only ♂spleen vs. ♂kidney.

In ♀liver vs. ♀spleen, ♀liver vs. ♀kidney and ♀spleen vs. ♀kidney, the q-values of 20, 16 and 6 pathways, respectively, were lower than 0.05. Only one category, namely, “glycine, serine and threonine metabolism”, was significantly enriched in ♀liver vs. ♀spleen, ♀liver vs. ♀kidney and ♀spleen vs. ♀kidney. “Biosynthesis of unsaturated fatty acids”, “carbon metabolism”, “cytokine-cytokine receptor interaction”, “fatty acid degradation”, “glycerolipid metabolism” and “peroxisome” were significantly enriched in only ♀liver vs. ♀spleen; “porphyrin and chlorophyll metabolism” and “renin-angiotensin system” were significantly enriched in only ♀liver vs. ♀kidney; “cytokine-cytokine receptor interaction”, “ECM-receptor interaction”, “focal adhesion”, “glutathione metabolism” and “primary bile acid biosynthesis” were significantly enriched in only ♀spleen vs. ♀kidney.

The enriched KEGG categories were also compared between male and female T. rubripes. In ♂liver vs. ♀liver, only two categories, namely, “PPAR signaling pathway” and “base excision repair”, were significantly enriched. In ♂spleen vs. ♀spleen, only one category, namely, “neuroactive ligand-receptor interaction”, was significantly enriched. In ♂kidney vs. ♀kidney, there was no KEGG category with a q-value < 0.05.

 

PPI network analysis

 

 

Fig. 4: Comparative Venn diagram of DEGs among different immune organs

 

 

 

Fig. 5: Scatter diagram of the top 20 enriched KEGG pathways of the DEGs from six groups. A red cross indicates that the q-value of the pathway is lower than 0.05

 

The PPI networks of DEGs were determined by BLAST analysis against the STRING database. In these PPI networks, there were 40, 50, 35, 44, 46, 37, 12, 15 and 4 modules, and 145,05, 20,167, 11,393, 17,492, 18,432, 11,955, 128, 113 and 66 nodes in liver vs. ♂spleen, ♂liver vs. ♂kidney, ♂spleen vs. ♂kidney, ♀liver vs. ♀spleen, ♀liver vs. ♀kidney, ♀spleen vs. ♀kidney, ♂liver vs. ♀liver, ♂spleen vs. ♀spleen and ♂kidney vs. ♀kidney, respectively.

 

Discussion

 

In this study, we performed a comparative RNA-Seq analysis of the important immune organs (livers, spleens and kidneys) of T. rubripes to find the key immune-related regulated networks, immune-related genes and DEGs of both the male and female fish using GO enrichment, KEGG pathway enrichment and PPI network analyses.

The most significantly enriched GO terms were “G-protein coupled receptor signaling pathway (GO:0007186)”, “integral component of membrane (GO:0016021)”, “G-protein coupled receptor activity (GO:0004930)” and “olfactory receptor activity (GO:0004984)” in the categories BP, CC and MF. Interestingly, all these GO terms are involved in cell surface receptor signaling. Cell surface receptors could bind cytokines, cell adhesion molecules and hormones and allow communication between cells and the external environment. In this study, a set of DEGs exhibited significant enrichment of these four GO terms. Previously, a transcriptomic study of common carp spleen infected with A. hydrophila also revealed that a number of genes associated with cell surface receptor signaling exhibit significantly different expression patterns after bacterial infection (Jiang et al. 2016). These results indicate that cell surface receptors are involved in immune-related signal transduction on cell surfaces and bind to extracellular ligands to activate, perpetuate, or inhibit an immune response.

In addition to GO enriched terms, we also found that the pathway “glycine, serine and threonine metabolism” was identified in all male and female groups (Fig. 5). Another important enriched pathway, namely, “cytokine-cytokine receptor interaction”, was identified in all the male groups, while in the female groups, this pathway was present in ♀liver vs. ♀spleen and ♀spleen vs. ♀kidney but not ♀liver vs. ♀kidney (Fig. 5). Cytokine-cytokine receptor interactions are involved in the innate immune response, a fundamental defense mechanism in fish. These interactions can stimulate the recruitment, activation, and adhesion of cells to sites of infection or injury (Laing and Secombes 2004). In this study, 107 DEGs, including chemokines and chemokine receptors (20 DEGs), ILs and IL receptors (15 DEGs), TNFs and TNF receptors (20 DEGs), were identified in this pathway (Table 3). Chemokines are a family of small cytokines and include the CC, CXC, CX3C and XC subfamilies. Cytokines exert their biological effects by interacting with chemokine receptors (Vicari et al. 1997). Members of the CXC chemokine, CC chemokine, CC receptor, and CXC receptor superfamilies have been identified from the channel catfish genome and are involved in defense against disease and hypoxia response (Fu et al. 2017a, b, c). Our previous study has also identified a set of chemokines and chemokine receptors in the gills of T. rubripes (Cui et al. 2014b). In tongue sole, 9 CC chemokines, including CCL20, CCL21, and CCL27, are cloned and characterized, and they act in the response to E. tarda and Vibrio harveyi (Hao and Li 2015). Four fish-specific CC chemokine receptors, namely, CCR4La, CCR4Lc1, CCR4Lc2 and CCR11, are identified in rainbow trout (Oncorhynchus mykiss), and the expression of these receptors change after challenges with Yersinia ruckeri (Qi et al. 2017). Meanwhile, CXC chemokines and CXC chemokine receptors are induced by the pathogens. CXC_F from Salmo trutta, osgCXCL12 and osgCXCR4 from Epinephelus coioides, and CXCR4 from Oplegnathus fasciatus are responsive to viral challenge with hemorrhagic septicemia virus and Y. ruckeri; nervous necrosis virus and E. tarda; S. iniae; and rock bream iridovirus, respectively (Thulasitha et al. 2015; Wu et al. 2015; Gorgoglione et al. 2016). ILs may trigger the protective defenses of the immune system and help eradicate pathogens (Jiang et al. 2016). In this study, 5 ILs and 10 IL receptor genes were identified (Table 3). A number of ILs and IL receptor genes have been identified in fish species. For example, IL-1β, IL-1RI and IL-1RII expression is induced by LPS and poly (I:C) in miiuy croaker, and overexpression of IL-1β increase IL-1RI expression (Yang et al. 2017). Seven IL17 ligands and four IL17 receptor homologs are identified and characterized from the transcriptomic and genomic databases for channel catfish, and three IL17A/Fs and their putative receptors, IL17RA and IL17RC, are induced following experimental challenge with E. ictaluri and F. columnare (Wang et al. 2014). Meanwhile, rock bream RbIL-15 and RbIL-15Rα are highly induced in the kidneys and spleen after infection with E. tarda, S. iniae and red sea bream iridovirus (Bae et al. 2013). TNF is another large class of inflammatory cytokines involved in cytokine-cytokine receptor interactions. In this study, 12 TNFs and 8 TNF receptor genes were identified (Table 3). In T. rubripes, ten novel TNFSF genes are identified and highly expressed in livers. In addition, the expression levels of three TNFSF10 genes are seen to be increased in poly (I:C)-treated head kidney cells (Biswas et al. 2015). The expressions of TNF-α and TNF-N in T. rubripes are significantly elevated by both the probiotic bacterial stimulants tested compared with the expression observed in the unstimulated control (Biswas et al. 2013). In addition, a TNF13B gene is cloned from Takifugu obscurus, and this gene may be a factor for the enhancement of immunological efficacy in fish (Ai et al. 2011).

DEG analysis between male and female T. rubripes showed that there were 5 DEGs in ♂liver vs. ♀liver, ♂spleen vs. ♀spleen and ♂kidney vs. ♀kidney (Fig. 4).
Table 3: Identification of DEGs in Cytokine-cytokine receptor interaction

 

Gene

♂liver vs. ♂spleen

♂liver vs. ♂kidney

♂spleen vs. ♂kidney

♀liver vs. ♀kidney

♀liver vs. ♀spleen

♀spleen vs. ♀kidney

Chemokine and chemokine receptor

LOC101061975 (CC25)

 

 

+

 

 

 

LOC101077134 (CC13)

 

 

 

+

 

+

ccl19l1

+

 

 

 

 

 

LOC105417015 (CC21)

 

 

 

 

 

+

cxcl11

 

 

 

 

 

+

cxcl14

 

 

 

+

 

 

cxcl12

 

 

 

 

 

+

ccl19

 

 

+

 

 

 

cxcl13

 

+

 

 

 

 

ccr10

+

+

+

+

 

+

ccr9

+

+

 

+

+

 

LOC101061903 (CCR4)

+

+

 

+

+

 

LOC101066745 (CCR9)

+

+

 

+

 

 

LOC101071120 (CCR4)

 

+

 

+

 

 

LOC101077667 (XCR1)

 

+

 

+

+

 

LOC101077892 (XCR1)

+

+

 

+

+

 

cxcr4

 

+

+

+

+

+

cxcr5

+

+

+

+

 

+

LOC101071867 (CXCR2)

 

+

+

+

+

 

LOC101079665 (CXCR)

 

 

 

+

+

 

TNF and TNF receptor

tnfa

+

 

+

 

 

 

tgfb2

+

 

 

 

 

+

tgfb3

 

 

 

+

 

 

LOC101066955 (TNFL6)

+

+

 

+

+

 

LOC101061420 (TNFL10)

+

+

+

+

+

 

LOC101075168 (TNFL13B)

+

 

 

+

 

 

tnfsf12

 

 

+

 

 

 

tnfsf11

+

+

 

+

 

 

tnfsf12

 

 

 

 

 

+

tnfsf13b

+

 

 

+

 

 

tnfsf10

+

+

 

 

+

 

LOC101061564 (TNF11B)

+

+

 

+

+

 

tnfrsf19

+

+

 

+

+

+

tnfrsf21

+

+

 

+

+

 

ngfr

+

 

+

+

 

+

edar

+

 

 

+

 

+

tgfbr2

 

 

 

+

 

 

tnfrsf11a

 

+

+

 

+

+

tnfrsf13b

+

+

+

+

+

+

VEGF and VEGF receptor

LOC101070662 (VEGFA)

 

+

+

 

+

+

vegfa

+

 

+

+

 

+

vegfb

+

+

+

+

 

+

vegfc

+

+

+

+

+

+

LOC101071220 (VEGFR)

+

 

+

+

 

+

flt1

 

 

 

 

 

+

flt3

+

+

 

+

+

+

flt4

 

 

 

 

 

+

kdr

 

 

 

 

 

+

BMP and BMP receptor

LOC101076733 (BMP7)

+

+

 

+

+

 

LOC101076832 (BMP7)

+

 

+

+

 

+

bmp7

+

+

 

+

+

 

bmpr1b

 

 

+

+

 

+

bmpr2

 

+

 

 

 

 

LOC101061264 (BMPR1B)

+

+

 

+

+

 

IL and IL receptor

il1b

 

 

 

+

 

 

il-18

+

 

+

+

 

+

il8

+

+

 

+

+

 

newGene_1959 (IL18)

+

 

+

+

 

+

il6st

 

+

 

+

+

 

il18rap

+

+

 

+

+

 

il22ra2

 

+

+

 

+

+

Table Continue


These 5 DEGs encoded complement C1q TNF-related protein 4-like, hepcidin-like, voltage-dependent calcium channel subunit alpha-2/delta-2-like, protein L-Myc-1b-like and apolipoprotein A-IV1 precursor. Table 3: Continued

 

il2rg

+

+

 

+

+

 

il7r

+

+

 

+

+

 

LOC101078915 (IL21)

+

+

 

+

+

 

LOC105417588 (IL31Rα)

 

+

+

+

+

+

LOC105419406 (IL2Rβ)

 

 

 

+

 

 

LOC778023 (IL8R1)

 

+

+

+

+

 

LOC101065602 (IL17RA)

+

+

 

+

+

 

LOC101075482 (IL12Rβ2)

+

 

 

 

 

+

Acv receptor

acvr1

+

 

+

 

 

 

acvr2b

 

 

 

+

 

 

LOC101064200 (acvr1b)

+

 

+

 

 

+

class I helical cytokine receptor

newGene_449 (class I CR)

+

+

+

+

+

+

newGene_1895 (class I CR)

+

+

+

+

+

+

LOC101075512 (class I CR)

+

+

 

+

+

 

PDGF

LOC101063539 (PDGFA)

 

 

 

+

+

 

pdgfa

 

 

+

 

 

 

pdgfra

 

 

+

 

 

 

LOC101074724 (PDGFRβ)

 

 

+

+

 

+

Others

LOC101061818 (LIFR)

+

+

 

+

+

 

LOC101065739 (AMCF2)

+

+

 

+

+

 

LOC101071352 (EPOR)

 

 

+

 

 

+

LOC101073748 (MRTPK)

+

+

 

+

+

+

LOC101079667 (PR-like)

+

 

+

+

+

+

LOC101069994 (M-CSF1R1)

+

+

+

+

+

 

LOC101080143 (M-CSF1R2)

 

 

 

 

+

 

LOC101075211 (G-CSF)

 

 

+

 

 

+

csf3r

 

+

+

+

+

 

LOC101077955 (TGFβ3)

+

 

 

+

 

 

LOC101078862 (HGF)

+

+

 

 

 

 

hgf

+

+

+

+

+

+

newGene_2384 (TPO)

+

+

 

+

+

 

mpl

+

 

+

+

 

+

egf

+

+

 

 

 

 

egfr

 

+

+

 

+

+

ghr1

+

+

 

+

+

 

kit

 

 

+

+

 

+

lepr

 

 

+

 

 

+

met

+

+

+

+

+

+

prlr

 

+

+

 

+

+

scf

+

 

 

 

 

 

LOC101067230

 

 

 

+

+

 

LOC101067572

+

+

+

+

+

+

LOC105417838

+

+

 

+

+

 

LOC105418164

+

+

 

+

+

 

LOC101075732

 

+

+

 

+

+

LOC101073224

 

+

 

+

+

 

 

Expression analysis showed that the expression of complement C1q TNF-related protein 4-like (CTRP4-like) was upregulated in all the male samples and the expression of voltage-dependent calcium channel subunit alpha-2/delta-2-like (CACNA2D2-like) and L-Myc-1b-like was upregulated in all the female samples. Complement C1q TNF-related protein 4 (CTRP4) is a secreted cytokine that is homologous to adiponectin, which plays an important role in immunity and metabolism (Huang et al. 2016). CTRP4 stimulates IL-6 synthesis and STAT3 and NFkB activation (Schäffler and Buechler 2012). Voltage-dependent calcium channels (VDCCs), also known as voltage-gated calcium channels (VGCCs), are a group of voltage-gated ion channels found in the membranes of excitable cells. In humans, it is found that VDCCs can mediate resistance to chemotherapy in small cell lung cancer (Yu et al. 2018). Human L-Myc is a proto-oncogenic transcription factor and acts by targeting the CCL6 chemokine gene. Mycl1is selectively expressed in the dendritic cells of the immune system and is controlled by IRF8 (Wumeshet al. 2014). An MYCL1 SNP, rs3134613, is associated with susceptibility to diffuse-type gastric cancer and with differentiation of gastric cancer in a southeast Chinese population (Chen et al. 2010).

 

Conclusion

 

In this study, RNA-Seq analysis of the livers, spleens and kidneys of T. rubripes was performed, and a number of DEGs were identified. The GO terms associated with cell surface receptor signaling were the most significantly enriched, and an immune-related pathway (cytokine-cytokine receptor interaction) was present in all the male groups. In addition, the DEGs between male and female T. rubripes were also analyzed, and CTRP4-like was found to be upregulated in all the male samples, while the expression of CACNA2D2-like and L-Myc-1b-like was upregulated in all the female samples. These results not only provide information regarding the molecular immune mechanisms in T. rubripes livers, spleens and kidneys but also identify candidates for breeding to enhance pathogen resistance in T. rubripes.

 

Acknowledgements

 

This study was supported by the National Key R&D Program of China (2018YFD0900301), the Science and Technology Innovation Fund Project of Dalian city (2018J12SN070), the High Level Talents Innovation Support Program of Dalian (2017RQ015) and the China Agriculture Research System (CARS-47).

 

References

 

Ai H, Y Shen, C Min, S Pang, J Zhang, S Zhang, Z Zhao (2011). Molecular structure, expression and bioactivity characterization of TNF13B (BAFF) gene in mefugu, Takifugu obscurus. Fish Shellf Immunol 30:1265‒1274

Asai K, N Hiki, Y Mimura, T Ogawa, K Unou, M Kaminishi (2001). Gender differences in cytokine secretion by human peripheral blood mononuclear cells: Role of estrogen in modulating LPS-induced cytokine secretion in an ex vivo septic model. Shock 16:340‒343

Bae JS, SH Shim, SD Hwang, JW Kim, DW Park, CI Park (2013). Molecular cloning and expression analysis of interleukin (IL)-15 and IL-15 receptor α from rock bream, Oplegnathus fasciatus. Fish Shellf Immunol 35:1209‒1215

Biswas G, S Kinoshita, T Kono, J Hikima, M Sakai (2015). Evolutionary evidence of tumor necrosis factor super family members in the Japanese pufferfish (Takifugu rubripes): Comprehensive genomic identification and expression analysis. Mar Genomics 22:25‒36

Biswas G, H Korenaga, R Nagamine, H Takayama, S Kawahara, S Takeda, Y Kikuchi, B Dashnyam, T Kono, M Sakai (2013). Cytokine responses in the Japanese pufferfish (Takifugu rubripes) head kidney cells induced with heat-killed probiotics isolated from the Mongolian dairy products. Fish Shellf Immunol 34:1170‒1177

Chen L, Z Liu, Y Su, D Wang, B Yin, B Shu, J Zhang, X Zhu, C Jia (2017). Characterization of Mycobacterium marinum infections in zebrafish wounds and sinus tracts. Wound Repair Regen 25:536‒540

Chen SQ, XD Lin, JW Zhu, Y Tang, JY Lin (2010). Association of a MYCL1 single nucleotide polymorphism, rs3134613, with susceptibility to diffuse-type gastric cancer and with differentiation of gastric cancer in a southeast Chinese population. DNA Cell Biol 29:739743

Cui J, N Jiang, X Hou, S Wu, Q Zhang, J Meng, Y Luan (2020). Genome-wide identification of lncRNAs and analysis of ceRNA networks during tomato resistance to Phytophthora infestans. Phytopathology 110:456‒464

Cui J, Y Luan, N Jiang, H Bao, J Meng (2017). Comparative transcriptome analysis between resistant and susceptible tomato allows the identification of lncRNA16397 conferring resistance to Phytophthora infestans by co-expressing glutaredoxin. Plant J 89:577‒589

Cui J, H Wang, S Liu, L Zhu, X Qiu, Z Jiang, X Wang, Z Liu (2014a). SNP discovery from transcriptome of the swim bladder of Takifugu rubripes. PLoS One 9; Article e92502

Cui J, S Liu, B Zhang, H Wang, H Sun, S Song, X Qiu, Y Liu, X Wang, Z Jiang, Z Liu (2014b). Transciptome analysis of the gill and swim bladder of Takifugu rubripes by RNA-Seq. PLoS One 9; Article e85505

Cui J, H Wang, S Liu, X Qiu, Z Jiang, X Wang (2014c). Transcriptome analysis of the gill of Takifugu rubripes using Illumina sequencing for discovery of SNPs. Compar Biochem Physiol D Genomics Proteomics 10:44‒51

Fu Q, Q Zeng, Y Li, Y Yang, C Li, S Liu, T Zhou, N Li, J Yao, C Jiang, D Li, Z Liu (2017a). The chemokinome superfamily in channel catfish: I. CXC subfamily and their involvement in disease defense and hypoxia responses. Fish Shellf Immunol 60:380‒390

Fu Q, Y Yang, C Li, Q Zeng, T Zhou, N Li, Y Liu, Y Li, X Wang, S Liu, D Li, Z Liu (2017b). The chemokinome superfamily: II. The 64 CC chemokines in channel catfish and their involvement in disease and hypoxia responses. Dev Compar Immunol 73:97‒108

Fu Q, Y Yang, C Li, Q Zeng, T Zhou, N Li, Y Liu, S Liu, Z Liu (2017c). The CC and CXC chemokine receptors in channel catfish (Ictalurus punctatus) and their involvement in disease and hypoxia responses. Dev Compar Immunol 77:241251

Gong Y, S Feng, S Li, Y Zhang, Z Zhao, M Hu, P Xu, Y Jiang (2017). Genome-wide characterization of Toll-like receptor gene family in common carp (Cyprinus carpio) and their involvement in host immune response to Aeromonas hydrophila infection. Compar Biochem Physiol D Genomics Proteomics 24:89‒98

Gorgoglione B, E Zahran, NG Taylor, SW Feist, J Zou, CJ Secombes (2016) Comparative study of CXC chemokines modulation in brown trout (Salmo trutta) following infection with a bacterial or viral pathogen. Mol Immunol 71:6477

Guo R, X Wang, H Su, X Zhang, H Liu (2017). Analysis and comparison of nutritional compositions in muscle, liver and skin of Takifugu rubripes. J Agric Univ Hebei 40:77‒82 (in Chinese)

Hamasaki M, Y Takeuchi, R Yazawa, S Yoshikawa, K Kadomura, T Yamada, K Miyaki, K Kikuchi, G Yoshizaki (2017). Production of tiger puffer Takifugu rubripes offspring from triploid grass puffer Takifugu niphobles parents. Mar Biotechnol 19:579‒591

Hao LX, MF Li (2015). Molecular characterization and expression analysis of nine CC chemokines in half-smooth tongue sole, Cynoglossus semilaevis. Fish Shellf Immunol 47:717724

Huang H, X Wu, L Cao, L Wang (2016). Preparation and Identification of Monoclonal Antibody Against C1q/TNF-Related Protein 4. Monoclon Antib Immunodiagn Immunother 35:280284

Jiang N, J Meng, J Cui, G Sun, Y Luan (2018). Function identification of miR482b, a negative regulator during tomato resistance to Phytophthora infestans. Hortic Res 5; Article 9

Jiang Y, S Feng, S Zhang, H Liu, J Feng, X Mu, X Sun, P Xu (2016). Transcriptome signatures in common carp spleen in response to Aeromonas hydrophila infection. Fish Shellf Immunol 57:41‒48

Kahlke V, C Dohm, K Brötzmann, S Schreiber, J Schröder (2000). Gender-related therapy: Early IL-10 administration after hemorrhage restores immune function in males but not in females. Shock 14:354‒359

Kanehisa M, M Araki, S Goto, M Hattori, M Hirakawa, M Itoh, T Katayama, S Kawashima, S Okuda, T Tokimatsu, Y Yamanishi (2008). KEGG for linking genomes to life and the environment. Nucl Acids Res 36:480‒484

Klein SL, I Marriott, EN Fish (2015). Sex-based differences in immune function and responses to vaccination. Trans Roy Soc Trop Med Hyg 109:9‒15

Kong D, Z Wang, L Zang, J Cui, H Li, H Sun, X Qiu, H Liu, X Wang (2019). Expression, purification and characterization of recombinant interferon-γ of Takifugu rubripes in Pichia pastoris. Intl J Agric Biol 21:19‒24

Kong D, J Cui, Z Wang, L Zang, H Sun, Z Hu, X Li, X Qiu, C Jiang, H Liu, T Zhang, S Liu, Z Jiang, X Meng, X Wang (2018). The regulatory networks conferred by IFN-γ in the kidney of Takifugu rubripes. Intl J Agric Biol 20:2189‒2195

Laing KJ, CJ Secombes (2004). Trout CC chemokines: Comparison of their sequences and expression patterns. Mol Immunol 41:793808

Ma A, W Li, X Wang, L Yue, Z Zhuang, X Meng, S Liu, L Tang, S Hou (2014). Research progress and outlook of Takifugu rubripes culture techniques. Mar Sci 38:116‒121

Mao X, T Cai, JG Olyarchuk, L Wei (2005). Automated genome annotation and pathway identification using the KEGG orthology (KO) as a controlled vocabulary. Bioinformatics 21:3787‒3793

Marriott I, KL Bost, YM Huet-Hudson (2006). Sexual dimorphism in expression of receptors for bacterial lipopolysaccharides in murine macrophages: A possible mechanism for gender-based differences in endotoxic shock susceptibility. J Reprod Immunol 71:12‒27

Moxley G, D Posthuma, P Carlson, E Estrada, J Han, LL Benson, MC Neale (2002). Sexual dimorphism in innate immunity. Arthritis Rheum 46:250‒258

Mu Y, F Ding, P Cui, J Ao, S Hu, X Chen (2010). Transcriptome and expression profiling analysis revealed changes of multiple signaling pathways involved in immunity in the large yellow croaker during Aeromonas hydrophila infection. BMC Genomics 11; Article 506

Ordas A, Z Hegedus, CV Henkel, OW Stockhammer, D Butler, HJ Jansen, P Racz, MMink, HP Spaink, AH Meijer (2011). Deep sequencing of the innate immune transcriptomic response of zebrafish embryos to Salmonella infection. Fish Shellf Immunol 31:716‒724

Peatman E, C Li, BC Peterson, DL Straus, BD Farmer, BH Beck (2013). Basal polarization of the mucosal compartment in Flavobacterium columnare susceptible and resistant channel catfish (Ictalurus punctatus). Mol Immunol 56:317‒327

Pennell LM, CL Galligan, EN Fish (2012). Sex affects immunity. J Autoimmun 38:282‒291

Qi Z, JW Holland, Y Jiang, CJ Secombes, P Nie, T Wang (2017). Molecular characterization and expression analysis of four fish-specific CC chemokine receptors CCR4La, CCR4Lc1, CCR4Lc2 and CCR11 in rainbow trout (Oncorhynchus mykiss). Fish Shellf Immunol 68:411427

Schäffler A, C Buechler (2012). CTRP family: Linking immunity to metabolism. Trends Endocrinol Metab 23:194204

Shannon P, A Markiel, O Ozier, NS Baliga, JT Wang, D Ramage, N Amin, B Schwikowski, T Ideker (2003). Cytoscape: A software environment for integrated models of biomolecular interaction networks. Genome Res 13:2498‒2504

Song X, X Hu, B Sun, Y Bo, K Wu, L Xiao, C Gong (2017). A transcriptome analysis focusing on inflammation-related genes of grass carp intestines following infection with Aeromonas hydrophila. Sci Rep 7; Article 40777

Stockhammer OW, H Rauwerda, FR Wittink, TM Breit, AH Meijer, HP Spaink (2010). Transcriptome analysis of Traf6 function in the innate immune response of zebrafish embryos. Mol Immunol 48:179‒190

Sudhagar A, G Kumar, M El-Matbouli (2018). Transcriptome analysis based on RNA-Seq in understanding pathogenic mechanisms of diseases and the immune system of fish: A comprehensive review. Intl J Mol Sci 19:245–271

Thulasitha WS, N Umasuthan, KS Revathy, I Whang, J Lee (2015). Molecular characterization, genomic structure and expressional profiles of a CXC chemokine receptor 4 (CXCR4) from rock bream Oplegnathus fasciatus. Fish Shellf Immunol 44:471477

Tran NT, ZX Gao, HHZhao, SK Yi, BX Chen, YH Zhao, L Lin, XQ Liu, WM Wang (2015). Transcriptome analysis and microsatellite discovery in the blunt snout bream (Megalobrama amblycephala) after challenge with Aeromonas hydrophila. Fish Shellf Immunol 45:72‒82

Trapnell C, L Pachter, SL Salzberg (2009). TopHat: Discovering splice junctions with RNA-Seq. Bioinformatics 25:1105‒1111

Vicari AP, DJ Figueroa, JA Hedrick, JS Foster, KP Singh, S Menon, NG Copeland, DJ Gilbert, NA Jenkins, KB Bacon, A Zlotnik (1997). TECK: A novel CC chemokine specifically expressed by thymic dendritic cells and potentially involved in T cell development. Immunity 7:291301

Wang B, Z Gan, S Cai, Z Wang, D Yu, Z Lin, Y Lu, Z Wu, J Jian (2016). Comprehensive identification and profiling of Nile tilapia (Oreochromis niloticus) microRNAs response to Streptococcus agalactiae infection through high-throughput sequencing. Fish Shellf Immunol 54:93‒106

Wang X, C Li, W Thongda, Y Luo, B Beck, E Peatman (2014). Characterization and mucosal responses of interleukin 17 family ligand and receptor genes in channel catfish Ictalurus punctatus. Fish Shellf Immunol 38:47‒55

Wang Z, J Cui, J Song, H Wang, K Gao, X Qiu, M Gou, X Li, Z Hu, X Wang, Y Chang (2018). Comparative Transcriptome analysis reveals growth-related genes in juvenile Chinese sea cucumber, Russian sea cucumber, and their hybrids. Mar Biotechnol 20:193‒205

Wu CS, TY Wang, CF Liu, HP Lin, YM Chen, TY Chen (2015). Molecular cloning and characterization of orange-spotted grouper (Epinephelus coioides) CXC. chemokine ligand 12. Fish Shellf Immunol 47:9961005

Wumesh KC, AT Satpathy, AS Rapaport, CG Briseño, X Wu, JC Albring, EV Russler-Germain, NM Kretzer, V Durai, SP Persaud, BT Edelson, J Loschko, M Cella, PM Allen, MC Nussenzweig, M Colonna, BP Sleckman, TL Murphy, KM Murphy (2014). L-Myc expression by dendritic cells is required for optimal T-cell priming. Nature 507:243247

Yang D, Q Liu, M Yang, H Wu, Q Wang, J Xiao, Y Zhang (2012). RNA-seq liver transcriptome analysis reveals an activated MHC-I pathway and an inhibited MHC-II pathway at the early stage of vaccine immunization in zebrafish. BMC Genomics 13; Article 319

Yang Q, Q Chu, X Zhao, T Xu (2017). Characterization of IL-1β and two types of IL-1 receptors in miiuy croaker and evolution analysis of IL-1 family. Fish Shellf Immunol 63:165‒172

Yang Y, H Yu, H Li, A Wang (2016). Transcriptome profiling of grass carp (Ctenopharyngodon idellus) infected with Aeromonas hydrophila. Fish Shellf Immunol 51:329‒336

Young MD, MJ Wakefield, GK Smyth, A Oshlack (2010). Gene ontology analysis for RNA-seq: Accounting for selection bias. Genome Biol 11:14–25

Yu J, S Wang, W Zhao, J Duan, Z Wang, H Chen, Y Tian, D Wang, J Zhao, T An, H Bai, M Wu, J Wang (2018). Mechanistic exploration of cancer stem cell marker voltage-dependent calcium channel α2δ1 subunit-mediated chemotherapy resistance in small cell lung cancer. Clin Cancer Res 24:21482158

Zhang X, S Wang, S Chen, Y Chen, Y Liu, C Shao, Q Wang, Y Lu, G Gong, S Ding, Z Sha (2015). Transcriptome analysis revealed changes of multiple genes involved in immunity in Cynoglossus semilaevis during Vibrio anguillarum infection. Fish Shellf Immunol 43:209‒218

Zhao C, M Fu, C Wang, Z Jiao, L Qiu (2016). RNA-Seq analysis of immune-relevant genes in Lateolabrax japonicus during Vibrio anguillarum infection. Fish Shellf Immunol 52:57‒64

Zhou X, J Cui, X Qiu, Y Chang, X Wang (2019). Whole genomic SNPs and SSRs development based on high-throughout transcript sequencing in sea cucumber (Apostichopus japonicus). Intl J Agric Biol 22:877‒881

Zhou X, J Cui, S Liu, D Kong, H Sun, C Gu, H Wang, X Qiu, Y Chang, Z Liu, X Wang (2016a). Comparative transcriptome analysis of papilla and skin in the sea cucumber, Apostichopus japonicus. Peer J 4; Article e1779

Zhou X, H Wang, J Cui, X Qiu, Y Chang, X Wang (2016b). Transcriptome analysis of tube foot and large scale marker discovery in sea cucumber, Apostichopus japonicus. Compar Biochem Physiol D Genomics Proteomics 20:4149

Zhu J, Q Fu, Q Ao, Y Tan, Y Luo, H Jiang, C Li, X Gan (2017). Transcriptomic profiling analysis of tilapia (Oreochromis niloticus) following Streptococcus agalactiae challenge. Fish Shellf Immunol 62:202‒212